# 1 "linpackd.f"
# 0 "linpackd.f"
C*$* Padding ( DD3, DD2, DD1 )
# 0 "linpackd.f"
C*$* Storage Order ( AA, A, B, DD1, X, DD2, TIME, DD3, IPVT )
*
*PLEASE NOTE THAT netlib HAS MOVED, THE NEW ADDRESS IS netlib@ornl.gov.
*THE OLD ADDRESS, netlib@mcs.anl.gov, WILL BE TURNED OFF SOON.
*
*** from netlib, Fri Jul 27 14:07:10 EDT 1990 ***
# 6 "linpackd.f"
      program MAIN
# 6 "linpackd.f"
      double precision second
      DOUBLEPRECISION AA(200,200), A(201,200), B(200), X(200)
# 8 "linpackd.f"
      DOUBLEPRECISION TIME(8,6), CRAY, OPS, TOTAL, NORMA, NORMX
# 9 "linpackd.f"
      double precision resid,residn,eps,epslon
      INTEGER IPVT(200)
# 11 "linpackd.f"
      INTEGER II1, II2, II3
# 11 "linpackd.f"
# 11 "linpackd.f"
# 11 "linpackd.f"
      DOUBLE PRECISION DD1 (129), DD2 (129), DD3 (665)
      LDA = 201
# 12 "linpackd.f"
      LDAA = 200
c
# 14 "linpackd.f"
      N = 100
# 15 "linpackd.f"
      cray = .056
      WRITE (6, 1)
# 17 "linpackd.f"
    1 format(' Please send the results of this run to:'//
     $       ' Jack J. Dongarra'/
     $       ' Computer Science Department'/
     $       ' University of Tennessee'/
     $       ' Knoxville, Tennessee 37996-1300'//
     $       ' Fax: 615-974-8296'//
     $       ' Internet: dongarra@cs.utk.edu'/)
      OPS = (2.D0 * 1000000) / 3.D0 + 2.D0 * 10000
c
# 26 "linpackd.f"
         call matgen(a,lda,n,b,norma)
         t1 = second()
         call dgefa(a,lda,n,ipvt,info)
         TIME(1,1) = SECOND(  ) - T1
# 30 "linpackd.f"
         t1 = second()
         CALL DGESL (A,LDA,N,IPVT,B,0)
# 32 "linpackd.f"
         TIME(1,2) = SECOND(  ) - T1
# 33 "linpackd.f"
         TOTAL = TIME(1,1) + TIME(1,2)
# 37 "linpackd.f"
         II1 = MOD (N, 4)
# 37 "linpackd.f"
         DO 3 I=1,II1
# 38 "linpackd.f"
            X(I) = B(I)
# 39 "linpackd.f"
    3       CONTINUE
c
c     compute a residual to verify results.
c
# 37 "linpackd.f"
C$DOACROSS IF(N .GT. 166),SHARE(II1,N,X,B),LOCAL(I)
# 37 "linpackd.f"
         DO 4 I=II1+1,N,4
# 38 "linpackd.f"
            x(i) = b(i)
            X(I+1) = B(I+1)
# 38 "linpackd.f"
            X(I+2) = B(I+2)
# 38 "linpackd.f"
            X(I+3) = B(I+3)
# 39 "linpackd.f"
    4    CONTINUE
# 40 "linpackd.f"
         call matgen(a,lda,n,b,norma)
# 41 "linpackd.f"
         II2 = MOD (N, 4)
# 41 "linpackd.f"
         DO 5 I=1,II2
# 42 "linpackd.f"
            B(I) = -B(I)
# 43 "linpackd.f"
    5       CONTINUE
# 41 "linpackd.f"
C$DOACROSS IF(N .GT. 125),SHARE(II2,N,B),LOCAL(I)
# 41 "linpackd.f"
         DO 6 I=II2+1,N,4
# 42 "linpackd.f"
            B(I) = -B(I)
# 42 "linpackd.f"
            B(I+1) = -B(I+1)
# 42 "linpackd.f"
            B(I+2) = -B(I+2)
# 42 "linpackd.f"
            B(I+3) = -B(I+3)
# 43 "linpackd.f"
    6    CONTINUE
# 44 "linpackd.f"
         call dmxpy(n,b,n,lda,x,a)
         RESID = 0.
# 46 "linpackd.f"
         NORMX = 0.
# 47 "linpackd.f"
         II3 = MOD (N, 4)
# 47 "linpackd.f"
         DO 2 I=1,II3
# 48 "linpackd.f"
            RESID = DMAX1 (RESID, DABS (B(I)))
# 49 "linpackd.f"
            NORMX = DMAX1 (NORMX, DABS (X(I)))
# 50 "linpackd.f"
    2    CONTINUE
# 47 "linpackd.f"
         DO 30 I=II3+1,N,4
# 48 "linpackd.f"
            resid = dmax1( resid, dabs(b(i)) )
            normx = dmax1( normx, dabs(x(i)) )
# 48 "linpackd.f"
            RESID = DMAX1 (RESID, DABS (B(I+1)))
# 49 "linpackd.f"
            NORMX = DMAX1 (NORMX, DABS (X(I+1)))
# 48 "linpackd.f"
            RESID = DMAX1 (RESID, DABS (B(I+2)))
# 49 "linpackd.f"
            NORMX = DMAX1 (NORMX, DABS (X(I+2)))
# 48 "linpackd.f"
            RESID = DMAX1 (RESID, DABS (B(I+3)))
# 49 "linpackd.f"
            NORMX = DMAX1 (NORMX, DABS (X(I+3)))
# 50 "linpackd.f"
   30    continue
         EPS = EPSLON( 1.D0 )
# 52 "linpackd.f"
         residn = resid/( n*norma*normx*eps )
         WRITE (6, 40)
# 54 "linpackd.f"
   40    format('     norm. resid      resid           machep',
     $          '         x(1)          x(n)')
         WRITE (6, 50) RESIDN, RESID, EPS, X(1), X(N)
# 57 "linpackd.f"
   50    format(1p5e16.8)
c
         WRITE (6, 60) N
# 60 "linpackd.f"
   60    format(//'    times are reported for matrices of order ',i5)
         WRITE (6, 70)
# 62 "linpackd.f"
   70    format(6x,'dgefa',6x,'dgesl',6x,'total',5x,'mflops',7x,'unit',
     $         6x,'ratio')
c
         TIME(1,3) = TOTAL
# 66 "linpackd.f"
         TIME(1,4) = OPS / (1.0D6 * TOTAL)
# 67 "linpackd.f"
         TIME(1,5) = 2.D0 / TIME(1,4)
# 68 "linpackd.f"
         TIME(1,6) = TOTAL / CRAY
# 69 "linpackd.f"
         WRITE (6, 80) LDA
# 70 "linpackd.f"
   80    format(' times for array with leading dimension of',i4)
         WRITE (6, 110) (TIME(1,I), I=1,6)
c
# 73 "linpackd.f"
         call matgen(a,lda,n,b,norma)
         t1 = second()
         call dgefa(a,lda,n,ipvt,info)
         TIME(2,1) = SECOND(  ) - T1
# 77 "linpackd.f"
         t1 = second()
         CALL DGESL (A,LDA,N,IPVT,B,0)
# 79 "linpackd.f"
         TIME(2,2) = SECOND(  ) - T1
# 80 "linpackd.f"
         TOTAL = TIME(2,1) + TIME(2,2)
# 81 "linpackd.f"
         TIME(2,3) = TOTAL
# 82 "linpackd.f"
         TIME(2,4) = OPS / (1.0D6 * TOTAL)
# 83 "linpackd.f"
         TIME(2,5) = 2.D0 / TIME(2,4)
# 84 "linpackd.f"
         TIME(2,6) = TOTAL / CRAY
c
# 86 "linpackd.f"
         call matgen(a,lda,n,b,norma)
         t1 = second()
         call dgefa(a,lda,n,ipvt,info)
         TIME(3,1) = SECOND(  ) - T1
# 90 "linpackd.f"
         t1 = second()
         CALL DGESL (A,LDA,N,IPVT,B,0)
# 92 "linpackd.f"
         TIME(3,2) = SECOND(  ) - T1
# 93 "linpackd.f"
         TOTAL = TIME(3,1) + TIME(3,2)
# 94 "linpackd.f"
         TIME(3,3) = TOTAL
# 95 "linpackd.f"
         TIME(3,4) = OPS / (1.0D6 * TOTAL)
# 96 "linpackd.f"
         TIME(3,5) = 2.D0 / TIME(3,4)
# 97 "linpackd.f"
         TIME(3,6) = TOTAL / CRAY
c
# 100 "linpackd.f"
         TM2 = 0
# 101 "linpackd.f"
         t1 = second()
         DO 90 I=1,10
# 103 "linpackd.f"
            tm = second()
            call matgen(a,lda,n,b,norma)
            tm2 = tm2 + second() - tm
            call dgefa(a,lda,n,ipvt,info)
   90    continue
         TIME(4,1) = (SECOND(  ) - T1 - TM2) / 10
# 109 "linpackd.f"
         t1 = second()
         DO 100 I=1,10
# 111 "linpackd.f"
            CALL DGESL (A,LDA,N,IPVT,B,0)
# 112 "linpackd.f"
  100    continue
         TIME(4,2) = (SECOND(  ) - T1) / 10
# 114 "linpackd.f"
         TOTAL = TIME(4,1) + TIME(4,2)
# 115 "linpackd.f"
         TIME(4,3) = TOTAL
# 116 "linpackd.f"
         TIME(4,4) = OPS / (1.0D6 * TOTAL)
# 117 "linpackd.f"
         TIME(4,5) = 2.D0 / TIME(4,4)
# 118 "linpackd.f"
         TIME(4,6) = TOTAL / CRAY
c
# 120 "linpackd.f"
         WRITE (6, 110) (TIME(2,I), I=1,6)
# 121 "linpackd.f"
         WRITE (6, 110) (TIME(3,I), I=1,6)
# 122 "linpackd.f"
         WRITE (6, 110) (TIME(4,I), I=1,6)
# 123 "linpackd.f"
  110    format(6(1pe11.3))
c
         call matgen(aa,ldaa,n,b,norma)
         t1 = second()
         call dgefa(aa,ldaa,n,ipvt,info)
         TIME(5,1) = SECOND(  ) - T1
# 129 "linpackd.f"
         t1 = second()
         CALL DGESL (AA,LDAA,N,IPVT,B,0)
# 131 "linpackd.f"
         TIME(5,2) = SECOND(  ) - T1
# 132 "linpackd.f"
         TOTAL = TIME(5,1) + TIME(5,2)
# 133 "linpackd.f"
         TIME(5,3) = TOTAL
# 134 "linpackd.f"
         TIME(5,4) = OPS / (1.0D6 * TOTAL)
# 135 "linpackd.f"
         TIME(5,5) = 2.D0 / TIME(5,4)
# 136 "linpackd.f"
         TIME(5,6) = TOTAL / CRAY
c
# 138 "linpackd.f"
         call matgen(aa,ldaa,n,b,norma)
         t1 = second()
         call dgefa(aa,ldaa,n,ipvt,info)
         TIME(6,1) = SECOND(  ) - T1
# 142 "linpackd.f"
         t1 = second()
         CALL DGESL (AA,LDAA,N,IPVT,B,0)
# 144 "linpackd.f"
         TIME(6,2) = SECOND(  ) - T1
# 145 "linpackd.f"
         TOTAL = TIME(6,1) + TIME(6,2)
# 146 "linpackd.f"
         TIME(6,3) = TOTAL
# 147 "linpackd.f"
         TIME(6,4) = OPS / (1.0D6 * TOTAL)
# 148 "linpackd.f"
         TIME(6,5) = 2.D0 / TIME(6,4)
# 149 "linpackd.f"
         TIME(6,6) = TOTAL / CRAY
c
# 151 "linpackd.f"
         call matgen(aa,ldaa,n,b,norma)
         t1 = second()
         call dgefa(aa,ldaa,n,ipvt,info)
         TIME(7,1) = SECOND(  ) - T1
# 155 "linpackd.f"
         t1 = second()
         CALL DGESL (AA,LDAA,N,IPVT,B,0)
# 157 "linpackd.f"
         TIME(7,2) = SECOND(  ) - T1
# 158 "linpackd.f"
         TOTAL = TIME(7,1) + TIME(7,2)
# 159 "linpackd.f"
         TIME(7,3) = TOTAL
# 160 "linpackd.f"
         TIME(7,4) = OPS / (1.0D6 * TOTAL)
# 161 "linpackd.f"
         TIME(7,5) = 2.D0 / TIME(7,4)
# 162 "linpackd.f"
         TIME(7,6) = TOTAL / CRAY
c
# 165 "linpackd.f"
         TM2 = 0
# 166 "linpackd.f"
         t1 = second()
         DO 120 I=1,10
# 168 "linpackd.f"
            tm = second()
            call matgen(aa,ldaa,n,b,norma)
            tm2 = tm2 + second() - tm
            call dgefa(aa,ldaa,n,ipvt,info)
  120    continue
         TIME(8,1) = (SECOND(  ) - T1 - TM2) / 10
# 174 "linpackd.f"
         t1 = second()
         DO 130 I=1,10
# 176 "linpackd.f"
            CALL DGESL (AA,LDAA,N,IPVT,B,0)
# 177 "linpackd.f"
  130    continue
         TIME(8,2) = (SECOND(  ) - T1) / 10
# 179 "linpackd.f"
         TOTAL = TIME(8,1) + TIME(8,2)
# 180 "linpackd.f"
         TIME(8,3) = TOTAL
# 181 "linpackd.f"
         TIME(8,4) = OPS / (1.0D6 * TOTAL)
# 182 "linpackd.f"
         TIME(8,5) = 2.D0 / TIME(8,4)
# 183 "linpackd.f"
         TIME(8,6) = TOTAL / CRAY
c
# 185 "linpackd.f"
         WRITE (6, 140) LDAA
# 186 "linpackd.f"
  140    format(/' times for array with leading dimension of',i4)
         WRITE (6, 110) (TIME(5,I), I=1,6)
# 188 "linpackd.f"
         WRITE (6, 110) (TIME(6,I), I=1,6)
# 189 "linpackd.f"
         WRITE (6, 110) (TIME(7,I), I=1,6)
# 190 "linpackd.f"
         WRITE (6, 110) (TIME(8,I), I=1,6)
# 191 "linpackd.f"
      stop
      end
# 193 "linpackd.f"
      subroutine matgen(a,lda,n,b,norma)
      DOUBLEPRECISION A(LDA,1), B(1), NORMA
# 195 "linpackd.f"
      INTEGER II1, II2, II3, II4
# 195 "linpackd.f"
      DOUBLE PRECISION DD1
# 195 "linpackd.f"
# 195 "linpackd.f"
c
      INIT = 1325
# 197 "linpackd.f"
      NORMA = 0.
# 198 "linpackd.f"
      DO 4 II1=1,N
# 199 "linpackd.f"
      II2 = MOD (N, 4)
# 199 "linpackd.f"
      DO 3 I=1,II2
# 200 "linpackd.f"
            INIT = MOD (INIT * 3125, 65536)
# 201 "linpackd.f"
            A(I,II1) = (INIT - 32768.) / 16384.
# 202 "linpackd.f"
            NORMA = DMAX1 (DABS (A(I,II1)), NORMA)
# 203 "linpackd.f"
    3       CONTINUE
# 199 "linpackd.f"
         DO 4 I=II2+1,N,4
# 200 "linpackd.f"
            INIT = MOD (INIT * 3125, 65536)
# 201 "linpackd.f"
            A(I,II1) = (INIT - 32768.) / 16384.
# 202 "linpackd.f"
            NORMA = DMAX1 (DABS (A(I,II1)), NORMA)
# 200 "linpackd.f"
            INIT = MOD (INIT * 3125, 65536)
# 201 "linpackd.f"
            A(I+1,II1) = (INIT - 32768.) / 16384.
# 202 "linpackd.f"
            NORMA = DMAX1 (DABS (A(I+1,II1)), NORMA)
# 200 "linpackd.f"
            INIT = MOD (INIT * 3125, 65536)
# 201 "linpackd.f"
            A(I+2,II1) = (INIT - 32768.) / 16384.
# 202 "linpackd.f"
            NORMA = DMAX1 (DABS (A(I+2,II1)), NORMA)
# 200 "linpackd.f"
            INIT = MOD (INIT * 3125, 65536)
# 201 "linpackd.f"
            A(I+3,II1) = (INIT - 32768.) / 16384.
# 202 "linpackd.f"
            NORMA = DMAX1 (DABS (A(I+3,II1)), NORMA)
# 203 "linpackd.f"
    4    CONTINUE
# 198 "linpackd.f"
         II3 = MOD (N, 4)
# 198 "linpackd.f"
         DO 5 II1=1,II3
# 206 "linpackd.f"
          B(II1) = 0.
# 207 "linpackd.f"
    5     CONTINUE
# 198 "linpackd.f"
C$DOACROSS IF(N .GT. 250),SHARE(II3,N,B),LOCAL(II1)
# 198 "linpackd.f"
          DO 6 II1=II3+1,N,4
# 206 "linpackd.f"
          B(II1) = 0.
# 206 "linpackd.f"
          B(II1+1) = 0.
# 206 "linpackd.f"
          B(II1+2) = 0.
# 206 "linpackd.f"
          B(II1+3) = 0.
# 207 "linpackd.f"
    6     CONTINUE
# 209 "linpackd.f"
C$DOACROSS IF(N .GT. 9),SHARE(N,B,A),LOCAL(DD1,II4,I,J)
# 209 "linpackd.f"
         DO 7 I=1,N
# 209 "linpackd.f"
         DD1 = B(I)
# 208 "linpackd.f"
         II4 = MOD (N, 4)
# 208 "linpackd.f"
         DO 2 J=1,II4
# 210 "linpackd.f"
            DD1 = DD1 + A(I,J)
# 212 "linpackd.f"
    2 CONTINUE
# 208 "linpackd.f"
      DO 50 J=II4+1,N,4
# 210 "linpackd.f"
            DD1 = DD1 + A(I,J)
# 210 "linpackd.f"
            DD1 = DD1 + A(I,J+1)
# 210 "linpackd.f"
            DD1 = DD1 + A(I,J+2)
# 210 "linpackd.f"
            DD1 = DD1 + A(I,J+3)
# 212 "linpackd.f"
   50 continue
      B(I) = DD1
# 211 "linpackd.f"
    7    CONTINUE
# 213 "linpackd.f"
      return
      end
# 215 "linpackd.f"
      subroutine dgefa(a,lda,n,ipvt,info)
      INTEGER LDA, N, IPVT(1), INFO
# 217 "linpackd.f"
      DOUBLEPRECISION A(LDA,1)
c
c     dgefa factors a double precision matrix by gaussian elimination.
c
c     dgefa is usually called by dgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
c
c     on entry
c
c        a       double precision(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that dgesl or dgedi will divide by zero
c                     if called.  use  rcond  in dgeco for a reliable
c                     indication of singularity.
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,dscal,idamax
c
c     internal variables
c
# 264 "linpackd.f"
      double precision t
      integer idamax,j,k,kp1,l,nm1
# 266 "linpackd.f"
      INTEGER II1, II2, II3, II4
# 266 "linpackd.f"
c
c
c     gaussian elimination with partial pivoting
c
      INFO = 0
# 273 "linpackd.f"
      DO 60 K=1,(N-1)
c
c        find l = pivot index
c
# 278 "linpackd.f"
         L = IDAMAX( N - K + 1, A(K,K), 1 ) + K - 1
# 279 "linpackd.f"
         ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
         IF (A(L,K) .NE. 0.D0) THEN
c
c           interchange if necessary
c
# 287 "linpackd.f"
            IF (L .NE. K) THEN
# 288 "linpackd.f"
               t = a(l,k)
               a(l,k) = a(k,k)
               a(k,k) = t
# 287 "linpackd.f"
               END IF
c
c           compute multipliers
c
# 295 "linpackd.f"
            T = -1.D0 / A(K,K)
# 296 "linpackd.f"
            CALL DSCAL (N - K,T,A(K+1,K),1)
# 302 "linpackd.f"
               IF (L .NE. K) THEN
# 302 "linpackd.f"
               II1 = K + 1
# 302 "linpackd.f"
               II2 = N - K
c
c           row elimination with column indexing
c
# 300 "linpackd.f"
            DO 30 J=II1,N
# 301 "linpackd.f"
               t = a(l,j)
                  a(l,j) = a(k,j)
                  a(k,j) = t
# 306 "linpackd.f"
               CALL DAXPY (II2,T,A(II1,K),1,A(II1,J),1)
# 307 "linpackd.f"
   30       continue
# 302 "linpackd.f"
            ELSE
# 302 "linpackd.f"
            II3 = K + 1
# 302 "linpackd.f"
            II4 = N - K
# 302 "linpackd.f"
            DO 4 J=II3,N
# 301 "linpackd.f"
               T = A(L,J)
# 306 "linpackd.f"
               CALL DAXPY (II4,T,A(II3,K),1,A(II3,J),1)
# 302 "linpackd.f"
    4          CONTINUE
# 302 "linpackd.f"
               END IF
# 283 "linpackd.f"
               ELSE
# 310 "linpackd.f"
            info = k
# 283 "linpackd.f"
            END IF
# 312 "linpackd.f"
   60 continue
# 314 "linpackd.f"
      ipvt(n) = n
      IF (A(N,N) .EQ. 0.D0) INFO = N
# 316 "linpackd.f"
      return
      end
# 318 "linpackd.f"
      subroutine dgesl(a,lda,n,ipvt,b,job)
      INTEGER LDA, N, IPVT(1), JOB
# 320 "linpackd.f"
      DOUBLEPRECISION A(LDA,1), B(1)
c
c     dgesl solves the double precision system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by dgeco or dgefa.
c
c     on entry
c
c        a       double precision(lda, n)
c                the output from dgeco or dgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from dgeco or dgefa.
c
c        b       double precision(n)
c                the right hand side vector.
c
c        job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  trans(a)*x = b  where
c                            trans(a)  is the transpose.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if dgeco has set rcond .gt. 0.0
c        or dgefa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call dgeco(a,lda,n,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call dgesl(a,lda,n,ipvt,c(1,j),0)
c        10 continue
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,ddot
c
c     internal variables
c
# 378 "linpackd.f"
      double precision ddot,t
      integer k,kb,l,nm1
c
# 382 "linpackd.f"
      IF (JOB .EQ. 0) THEN
c
c        job = 0 , solve  a * x = b
c        first solve  l*y = b
c
# 388 "linpackd.f"
         DO 20 K=1,(N-1)
# 390 "linpackd.f"
            T = B(IPVT(K))
# 391 "linpackd.f"
            IF (IPVT(K) .NE. K) THEN
# 392 "linpackd.f"
               B(IPVT(K)) = B(K)
# 393 "linpackd.f"
               b(k) = t
# 391 "linpackd.f"
               END IF
# 395 "linpackd.f"
            CALL DAXPY (N - K,T,A(K+1,K),1,B(K+1),1)
# 396 "linpackd.f"
   20    continue
c
c        now solve  u*x = y
c
# 401 "linpackd.f"
         DO 40 KB=1,N
# 402 "linpackd.f"
            K = N - KB + 1
# 403 "linpackd.f"
            b(k) = b(k)/a(k,k)
            T = -B(K)
# 405 "linpackd.f"
            CALL DAXPY (K - 1,T,A(1,K),1,B(1),1)
# 406 "linpackd.f"
   40    continue
# 382 "linpackd.f"
         ELSE
c
c        job = nonzero, solve  trans(a) * x = b
c        first solve  trans(u)*y = b
c
# 413 "linpackd.f"
         DO 60 K=1,N
# 414 "linpackd.f"
            T = DDOT( K - 1, A(1,K), 1, B(1), 1 )
# 415 "linpackd.f"
            b(k) = (b(k) - t)/a(k,k)
   60    continue
c
c        now solve trans(l)*x = y
c
# 421 "linpackd.f"
         DO 80 KB=1,(N-1)
# 423 "linpackd.f"
            B((N-KB)) = B((N-KB)) + DDOT( N - (N - KB), A((N-KB)+1,(N-KB
     X        )), 1, B((N-KB)+1), 1 )
# 425 "linpackd.f"
            IF (IPVT((N-KB)) .NE. (N - KB)) THEN
# 426 "linpackd.f"
               T = B(IPVT((N-KB)))
# 427 "linpackd.f"
               B(IPVT((N-KB))) = B((N-KB))
# 428 "linpackd.f"
               B((N-KB)) = T
# 425 "linpackd.f"
               END IF
# 430 "linpackd.f"
   80    continue
# 382 "linpackd.f"
         END IF
# 433 "linpackd.f"
      return
      end
# 435 "linpackd.f"
      subroutine daxpy(n,da,dx,incx,dy,incy)
c
c     constant times a vector plus a vector.
c     jack dongarra, linpack, 3/11/78.
c
      DOUBLEPRECISION DX(1), DY(1), DA
# 441 "linpackd.f"
      integer i,incx,incy,ix,iy,m,mp1,n
# 442 "linpackd.f"
      INTEGER II1, II2, II3, II4, II5, II6
# 442 "linpackd.f"
c
      IF (N .LE. 0) RETURN
# 444 "linpackd.f"
      IF (DA .EQ. 0.D0) RETURN
# 445 "linpackd.f"
      IF (INCX .NE. 1 .OR. INCY .NE. 1) THEN
c
c        code for unequal increments or equal increments
c          not equal to 1
c
# 450 "linpackd.f"
      IX = 1
# 451 "linpackd.f"
      IY = 1
# 452 "linpackd.f"
      IF (INCX .LT. 0) IX = (1 - N) * INCX + 1
# 453 "linpackd.f"
      IF (INCY .LT. 0) IY = (1 - N) * INCY + 1
# 457 "linpackd.f"
      II1 = IY
# 456 "linpackd.f"
      II2 = IX
# 454 "linpackd.f"
      IF (INCY .EQ. 0) THEN
# 454 "linpackd.f"
      II4 = MOD (N, 4)
# 454 "linpackd.f"
      DO 3 I=1,II4
# 455 "linpackd.f"
        DY(IY) = DY(IY) + DA * DX(IX)
# 456 "linpackd.f"
        IX = IX + INCX
# 457 "linpackd.f"
        IY = IY + INCY
# 458 "linpackd.f"
    3 CONTINUE
# 454 "linpackd.f"
      DO 10 I=II4+1,N,4
# 455 "linpackd.f"
        dy(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
# 455 "linpackd.f"
        DY(IY) = DY(IY) + DA * DX(IX)
# 456 "linpackd.f"
        IX = IX + INCX
# 457 "linpackd.f"
        IY = IY + INCY
# 455 "linpackd.f"
        DY(IY) = DY(IY) + DA * DX(IX)
# 456 "linpackd.f"
        IX = IX + INCX
# 457 "linpackd.f"
        IY = IY + INCY
# 455 "linpackd.f"
        DY(IY) = DY(IY) + DA * DX(IX)
# 456 "linpackd.f"
        IX = IX + INCX
# 457 "linpackd.f"
        IY = IY + INCY
# 458 "linpackd.f"
   10 continue
# 454 "linpackd.f"
      ELSE
# 454 "linpackd.f"
      II5 = MOD (N, 4)
# 454 "linpackd.f"
      DO 4 I=1,II5
# 455 "linpackd.f"
        II3 = II1 + (I - 1) * INCY
# 455 "linpackd.f"
        DY(II3) = DY(II3) + DA * DX(II2+(I-1)*INCX)
# 458 "linpackd.f"
    4 CONTINUE
# 454 "linpackd.f"
C$DOACROSS IF(N .GT. 25),SHARE(II5,N,II1,INCY,DY,DA,II2,INCX,DX),LOCAL(
C$&  II3,I)
# 454 "linpackd.f"
      DO 2 I=II5+1,N,4
# 455 "linpackd.f"
        II3 = II1 + (I - 1) * INCY
# 455 "linpackd.f"
        DY(II3) = DY(II3) + DA * DX(II2+(I-1)*INCX)
# 455 "linpackd.f"
        II3 = II1 + I * INCY
# 455 "linpackd.f"
        DY(II3) = DY(II3) + DA * DX(II2+I*INCX)
# 455 "linpackd.f"
        II3 = II1 + (I + 1) * INCY
# 455 "linpackd.f"
        DY(II3) = DY(II3) + DA * DX(II2+(I+1)*INCX)
# 455 "linpackd.f"
        II3 = II1 + (I + 2) * INCY
# 455 "linpackd.f"
        DY(II3) = DY(II3) + DA * DX(II2+(I+2)*INCX)
# 458 "linpackd.f"
    2 CONTINUE
# 454 "linpackd.f"
      END IF
# 459 "linpackd.f"
      return
# 445 "linpackd.f"
      END IF
c
c        code for both increments equal to 1
c
# 464 "linpackd.f"
      II6 = MOD (N, 4)
# 464 "linpackd.f"
      DO 5 I=1,II6
# 465 "linpackd.f"
        DY(I) = DY(I) + DA * DX(I)
# 466 "linpackd.f"
    5   CONTINUE
# 464 "linpackd.f"
C$DOACROSS IF(N .GT. 83),SHARE(II6,N,DY,DA,DX),LOCAL(I)
# 464 "linpackd.f"
      DO 6 I=II6+1,N,4
# 465 "linpackd.f"
        dy(i) = dy(i) + da*dx(i)
        DY(I+1) = DY(I+1) + DA * DX(I+1)
# 465 "linpackd.f"
        DY(I+2) = DY(I+2) + DA * DX(I+2)
# 465 "linpackd.f"
        DY(I+3) = DY(I+3) + DA * DX(I+3)
# 466 "linpackd.f"
    6 CONTINUE
# 467 "linpackd.f"
      return
      end
# 469 "linpackd.f"
      double precision function ddot(n,dx,incx,dy,incy)
c
c     forms the dot product of two vectors.
c     jack dongarra, linpack, 3/11/78.
c
      DOUBLEPRECISION DX(1), DY(1), DTEMP
# 475 "linpackd.f"
      integer i,incx,incy,ix,iy,m,mp1,n
# 476 "linpackd.f"
      INTEGER II3, II4
# 476 "linpackd.f"
c
      DDOT = 0.D0
# 478 "linpackd.f"
      DTEMP = 0.D0
# 479 "linpackd.f"
      IF (N .LE. 0) RETURN
# 480 "linpackd.f"
      IF (INCX .NE. 1 .OR. INCY .NE. 1) THEN
c
c        code for unequal increments or equal increments
c          not equal to 1
c
# 485 "linpackd.f"
      IX = 1
# 486 "linpackd.f"
      IY = 1
# 487 "linpackd.f"
      IF (INCX .LT. 0) IX = (1 - N) * INCX + 1
# 488 "linpackd.f"
      IF (INCY .LT. 0) IY = (1 - N) * INCY + 1
# 489 "linpackd.f"
      II3 = MOD (N, 4)
# 489 "linpackd.f"
      DO 2 I=1,II3
# 490 "linpackd.f"
        DTEMP = DTEMP + DX(IX) * DY(IY)
# 491 "linpackd.f"
        IX = IX + INCX
# 492 "linpackd.f"
        IY = IY + INCY
# 493 "linpackd.f"
    2 CONTINUE
# 489 "linpackd.f"
      DO 10 I=II3+1,N,4
# 490 "linpackd.f"
        dtemp = dtemp + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
# 490 "linpackd.f"
        DTEMP = DTEMP + DX(IX) * DY(IY)
# 491 "linpackd.f"
        IX = IX + INCX
# 492 "linpackd.f"
        IY = IY + INCY
# 490 "linpackd.f"
        DTEMP = DTEMP + DX(IX) * DY(IY)
# 491 "linpackd.f"
        IX = IX + INCX
# 492 "linpackd.f"
        IY = IY + INCY
# 490 "linpackd.f"
        DTEMP = DTEMP + DX(IX) * DY(IY)
# 491 "linpackd.f"
        IX = IX + INCX
# 492 "linpackd.f"
        IY = IY + INCY
# 493 "linpackd.f"
   10 continue
      ddot = dtemp
      return
# 480 "linpackd.f"
      END IF
c
c        code for both increments equal to 1
c
# 500 "linpackd.f"
      II4 = MOD (N, 4)
# 500 "linpackd.f"
      DO 3 I=1,II4
# 501 "linpackd.f"
        DTEMP = DTEMP + DX(I) * DY(I)
# 502 "linpackd.f"
    3 CONTINUE
# 500 "linpackd.f"
      DO 30 I=II4+1,N,4
# 501 "linpackd.f"
        dtemp = dtemp + dx(i)*dy(i)
        DTEMP = DTEMP + DX(I+1) * DY(I+1)
# 501 "linpackd.f"
        DTEMP = DTEMP + DX(I+2) * DY(I+2)
# 501 "linpackd.f"
        DTEMP = DTEMP + DX(I+3) * DY(I+3)
# 502 "linpackd.f"
   30 continue
      ddot = dtemp
      return
      end
# 506 "linpackd.f"
      subroutine  dscal(n,da,dx,incx)
c
c     scales a vector by a constant.
c     jack dongarra, linpack, 3/11/78.
c
      DOUBLEPRECISION DA, DX(1)
# 512 "linpackd.f"
      integer i,incx,m,mp1,n,nincx
# 513 "linpackd.f"
      INTEGER II1, II2
# 513 "linpackd.f"
c
      IF (N .LE. 0) RETURN
# 515 "linpackd.f"
      IF (INCX .NE. 1) THEN
c
c        code for increment not equal to 1
c
# 520 "linpackd.f"
      II1 = N * INCX
# 520 "linpackd.f"
C$DOACROSS IF(II1 / INCX .GT. 125),SHARE(II1,INCX,DX,DA),LOCAL(I)
# 520 "linpackd.f"
      DO 2 I=1,II1,INCX
# 521 "linpackd.f"
        dx(i) = da*dx(i)
# 522 "linpackd.f"
    2 CONTINUE
# 523 "linpackd.f"
      return
# 515 "linpackd.f"
      END IF
c
c        code for increment equal to 1
c
# 528 "linpackd.f"
      II2 = MOD (N, 4)
# 528 "linpackd.f"
      DO 3 I=1,II2
# 529 "linpackd.f"
        DX(I) = DA * DX(I)
# 530 "linpackd.f"
    3   CONTINUE
# 528 "linpackd.f"
C$DOACROSS IF(N .GT. 125),SHARE(II2,N,DX,DA),LOCAL(I)
# 528 "linpackd.f"
      DO 4 I=II2+1,N,4
# 529 "linpackd.f"
        dx(i) = da*dx(i)
        DX(I+1) = DA * DX(I+1)
# 529 "linpackd.f"
        DX(I+2) = DA * DX(I+2)
# 529 "linpackd.f"
        DX(I+3) = DA * DX(I+3)
# 530 "linpackd.f"
    4 CONTINUE
# 531 "linpackd.f"
      return
      end
# 533 "linpackd.f"
      integer function idamax(n,dx,incx)
c
c     finds the index of element having max. dabsolute value.
c     jack dongarra, linpack, 3/11/78.
c
      DOUBLEPRECISION DX(1), DMAX
# 539 "linpackd.f"
      integer i,incx,ix,n
# 540 "linpackd.f"
      INTEGER II1, II2, II3
# 540 "linpackd.f"
c
      IDAMAX = 0
# 542 "linpackd.f"
      IF (N .LT. 1) RETURN
# 543 "linpackd.f"
      IDAMAX = 1
# 544 "linpackd.f"
      IF (N .EQ. 1) RETURN
# 545 "linpackd.f"
      IF (INCX .NE. 1) THEN
c
c        code for increment not equal to 1
c
# 550 "linpackd.f"
      DMAX = DABS (DX(1))
# 551 "linpackd.f"
      IX = INCX + 1
# 552 "linpackd.f"
      II2 = MOD (N - 1, 4)
# 552 "linpackd.f"
      DO 2 I=2,II2+1
# 553 "linpackd.f"
         IF (DABS (DX(IX)) .GT. DMAX) THEN
# 554 "linpackd.f"
         IDAMAX = I
# 555 "linpackd.f"
         DMAX = DABS (DX(IX))
# 553 "linpackd.f"
         END IF
# 556 "linpackd.f"
         IX = IX + INCX
# 557 "linpackd.f"
    2 CONTINUE
# 552 "linpackd.f"
      DO 10 I=II2+2,N,4
# 553 "linpackd.f"
         IF (DABS (DX(IX)) .GT. DMAX) THEN
# 554 "linpackd.f"
         idamax = i
         DMAX = DABS (DX(IX))
# 553 "linpackd.f"
         END IF
# 556 "linpackd.f"
         IX = IX + INCX
# 553 "linpackd.f"
         IF (DABS (DX(IX)) .GT. DMAX) THEN
# 554 "linpackd.f"
         IDAMAX = I + 1
# 555 "linpackd.f"
         DMAX = DABS (DX(IX))
# 553 "linpackd.f"
         END IF
# 556 "linpackd.f"
         IX = IX + INCX
# 553 "linpackd.f"
         IF (DABS (DX(IX)) .GT. DMAX) THEN
# 554 "linpackd.f"
         IDAMAX = I + 2
# 555 "linpackd.f"
         DMAX = DABS (DX(IX))
# 553 "linpackd.f"
         END IF
# 556 "linpackd.f"
         IX = IX + INCX
# 553 "linpackd.f"
         IF (DABS (DX(IX)) .GT. DMAX) THEN
# 554 "linpackd.f"
         IDAMAX = I + 3
# 555 "linpackd.f"
         DMAX = DABS (DX(IX))
# 553 "linpackd.f"
         END IF
# 556 "linpackd.f"
         IX = IX + INCX
# 557 "linpackd.f"
   10 continue
      return
# 545 "linpackd.f"
      END IF
c
c        code for increment equal to 1
c
# 562 "linpackd.f"
      DMAX = DABS (DX(1))
# 563 "linpackd.f"
      II3 = MOD (N - 1, 4)
# 563 "linpackd.f"
      DO 3 I=2,II3+1
# 564 "linpackd.f"
         IF (DABS (DX(I)) .GT. DMAX) THEN
# 565 "linpackd.f"
         IDAMAX = I
# 566 "linpackd.f"
         DMAX = DABS (DX(I))
# 564 "linpackd.f"
         END IF
# 567 "linpackd.f"
    3 CONTINUE
# 563 "linpackd.f"
      DO 30 I=II3+2,N,4
# 564 "linpackd.f"
         IF (DABS (DX(I)) .GT. DMAX) THEN
# 565 "linpackd.f"
         idamax = i
         dmax = dabs(dx(i))
# 564 "linpackd.f"
         END IF
# 564 "linpackd.f"
         IF (DABS (DX(I+1)) .GT. DMAX) THEN
# 565 "linpackd.f"
         IDAMAX = I + 1
# 566 "linpackd.f"
         DMAX = DABS (DX(I+1))
# 564 "linpackd.f"
         END IF
# 564 "linpackd.f"
         IF (DABS (DX(I+2)) .GT. DMAX) THEN
# 565 "linpackd.f"
         IDAMAX = I + 2
# 566 "linpackd.f"
         DMAX = DABS (DX(I+2))
# 564 "linpackd.f"
         END IF
# 564 "linpackd.f"
         IF (DABS (DX(I+3)) .GT. DMAX) THEN
# 565 "linpackd.f"
         IDAMAX = I + 3
# 566 "linpackd.f"
         DMAX = DABS (DX(I+3))
# 564 "linpackd.f"
         END IF
# 567 "linpackd.f"
   30 continue
      return
      end
# 570 "linpackd.f"
      double precision function epslon (x)
      double precision x
c
c     estimate unit roundoff in quantities of size x.
c
      double precision a,b,c,eps
# 576 "linpackd.f"
      INTEGER II1
# 576 "linpackd.f"
      DOUBLE PRECISION DD1
# 576 "linpackd.f"
# 576 "linpackd.f"
c
c     this program should function properly on all systems
c     satisfying the following two assumptions,
c        1.  the base used in representing dfloating point
c            numbers is not a power of three.
c        2.  the quantity  a  in statement 10 is represented to
c            the accuracy used in dfloating point variables
c            that are stored in memory.
c     the statement number 10 and the go to 10 are intended to
c     force optimizing compilers to generate code satisfying
c     assumption 2.
c     under these assumptions, it should be true that,
c            a  is not exactly equal to four-thirds,
c            b  has a zero for its last bit or digit,
c            c  is not exactly equal to one,
c            eps  measures the separation of 1.0 from
c                 the next larger dfloating point number.
c     the developers of eispack would appreciate being informed
c     about any systems where these assumptions do not hold.
c
c     *****************************************************************
c     this routine is one of the auxiliary routines used by eispack iii
c     to avoid machine dependencies.
c     *****************************************************************
c
c     this version dated 4/6/83.
c
      A = 4.D0 / 3.D0
# 604 "linpackd.f"
      II1 = 1
# 604 "linpackd.f"
      DD1 = A - 1.D0
# 604 "linpackd.f"
      DO   WHILE ( II1 .EQ. 1 )
# 604 "linpackd.f"
      B = DD1
# 605 "linpackd.f"
      c = b + b + b
      EPS = DABS (C - 1.D0)
# 607 "linpackd.f"
      IF (EPS .NE. 0.D0) II1 = 0
# 607 "linpackd.f"
      END DO
# 608 "linpackd.f"
      epslon = eps*dabs(x)
      return
      end
# 611 "linpackd.f"
      subroutine dmxpy (n1, y, n2, ldm, x, m)
      double precision y(*), x(*), m(ldm,*)
# 613 "linpackd.f"
      INTEGER II1, II2, II3, II4, II5, II6, II7, II8
# 613 "linpackd.f"
      DOUBLE PRECISION DD1, DD2, DD3, DD4, DD5, DD6, DD7, DD8, DD9, DD10
     X  , DD11, DD12, DD13, DD14, DD15, DD16
# 613 "linpackd.f"
# 613 "linpackd.f"
c
c   purpose:
c     multiply matrix m times vector x and add the result to vector y.
c
c   parameters:
c
c     n1 integer, number of elements in vector y, and number of rows in
c         matrix m
c
c     y double precision(n1), vector of length n1 to which is added
c         the product m*x
c
c     n2 integer, number of elements in vector x, and number of columns
c         in matrix m
c
c     ldm integer, leading dimension of array m
c
c     x double precision(n2), vector of length n2
c
c     m double precision(ldm,n2), matrix of n1 rows and n2 columns
c
c ----------------------------------------------------------------------
c
c   cleanup odd vector
c
      IF (MOD (N2, 2) .GE. 1) THEN
# 640 "linpackd.f"
      II1 = MOD (N2, 2)
# 640 "linpackd.f"
      DD1 = X(II1)
# 640 "linpackd.f"
      II5 = MOD (N1, 4)
# 640 "linpackd.f"
      DO 2 I=1,II5
# 641 "linpackd.f"
            Y(I) = Y(I) + DD1 * M(I,II1)
# 642 "linpackd.f"
    2       CONTINUE
# 640 "linpackd.f"
C$DOACROSS IF(N1 .GT. 41),SHARE(II5,N1,Y,DD1,II1,M),LOCAL(I)
# 640 "linpackd.f"
         DO 3 I=II5+1,N1,4
# 641 "linpackd.f"
            Y(I) = Y(I) + DD1 * M(I,II1)
# 641 "linpackd.f"
            Y(I+1) = Y(I+1) + DD1 * M(I+1,II1)
# 641 "linpackd.f"
            Y(I+2) = Y(I+2) + DD1 * M(I+2,II1)
# 641 "linpackd.f"
            Y(I+3) = Y(I+3) + DD1 * M(I+3,II1)
# 642 "linpackd.f"
    3    CONTINUE
# 643 "linpackd.f"
      endif
c
c   cleanup odd group of two vectors
c
# 648 "linpackd.f"
      IF (MOD (N2, 4) .GE. 2) THEN
# 649 "linpackd.f"
      II2 = MOD (N2, 4)
# 649 "linpackd.f"
      DD2 = X(II2)
# 649 "linpackd.f"
      DD3 = X(II2-1)
# 649 "linpackd.f"
      II6 = MOD (N1, 4)
# 649 "linpackd.f"
      DO 4 I=1,II6
# 650 "linpackd.f"
            Y(I) = (Y(I) + DD3 * M(I,II2-1)) + DD2 * M(I,II2)
# 652 "linpackd.f"
    4       CONTINUE
# 649 "linpackd.f"
C$DOACROSS IF(N1 .GT. 21),SHARE(II6,N1,Y,DD3,II2,M,DD2),LOCAL(I)
# 649 "linpackd.f"
         DO 5 I=II6+1,N1,4
# 650 "linpackd.f"
            Y(I) = (Y(I) + DD3 * M(I,II2-1)) + DD2 * M(I,II2)
# 650 "linpackd.f"
            Y(I+1) = (Y(I+1) + DD3 * M(I+1,II2-1)) + DD2 * M(I+1,II2)
# 650 "linpackd.f"
            Y(I+2) = (Y(I+2) + DD3 * M(I+2,II2-1)) + DD2 * M(I+2,II2)
# 650 "linpackd.f"
            Y(I+3) = (Y(I+3) + DD3 * M(I+3,II2-1)) + DD2 * M(I+3,II2)
# 652 "linpackd.f"
    5    CONTINUE
# 653 "linpackd.f"
      endif
c
c   cleanup odd group of four vectors
c
# 658 "linpackd.f"
      IF (MOD (N2, 8) .GE. 4) THEN
# 659 "linpackd.f"
      II3 = MOD (N2, 8)
# 659 "linpackd.f"
      DD4 = X(II3)
# 659 "linpackd.f"
      DD5 = X(II3-1)
# 659 "linpackd.f"
      DD6 = X(II3-2)
# 659 "linpackd.f"
      DD7 = X(II3-3)
# 659 "linpackd.f"
      II7 = MOD (N1, 4)
# 659 "linpackd.f"
      DO 6 I=1,II7
# 660 "linpackd.f"
            Y(I) = (((Y(I) + DD7 * M(I,II3-3)) + DD6 * M(I,II3-2)) + DD5
     X              * M(I,II3-1)) + DD4 * M(I,II3)
# 663 "linpackd.f"
    6       CONTINUE
# 659 "linpackd.f"
C$DOACROSS IF(N1 .GT. 11),SHARE(II7,N1,Y,DD7,II3,M,DD6,DD5,DD4),LOCAL(I)
# 659 "linpackd.f"
         DO 7 I=II7+1,N1,4
# 660 "linpackd.f"
            Y(I) = (((Y(I) + DD7 * M(I,II3-3)) + DD6 * M(I,II3-2)) + DD5
     X              * M(I,II3-1)) + DD4 * M(I,II3)
# 660 "linpackd.f"
            Y(I+1) = (((Y(I+1) + DD7 * M(I+1,II3-3)) + DD6 * M(I+1,II3-2
     X             )) + DD5 * M(I+1,II3-1)) + DD4 * M(I+1,II3)
# 660 "linpackd.f"
            Y(I+2) = (((Y(I+2) + DD7 * M(I+2,II3-3)) + DD6 * M(I+2,II3-2
     X             )) + DD5 * M(I+2,II3-1)) + DD4 * M(I+2,II3)
# 660 "linpackd.f"
            Y(I+3) = (((Y(I+3) + DD7 * M(I+3,II3-3)) + DD6 * M(I+3,II3-2
     X             )) + DD5 * M(I+3,II3-1)) + DD4 * M(I+3,II3)
# 663 "linpackd.f"
    7    CONTINUE
# 664 "linpackd.f"
      endif
c
c   cleanup odd group of eight vectors
c
# 669 "linpackd.f"
      IF (MOD (N2, 16) .GE. 8) THEN
# 670 "linpackd.f"
      II4 = MOD (N2, 16)
# 670 "linpackd.f"
      DD8 = X(II4)
# 670 "linpackd.f"
      DD9 = X(II4-1)
# 670 "linpackd.f"
      DD10 = X(II4-2)
# 670 "linpackd.f"
      DD11 = X(II4-3)
# 670 "linpackd.f"
      DD12 = X(II4-4)
# 670 "linpackd.f"
      DD13 = X(II4-5)
# 670 "linpackd.f"
      DD14 = X(II4-6)
# 670 "linpackd.f"
      DD15 = X(II4-7)
# 670 "linpackd.f"
      II8 = MOD (N1, 2)
# 670 "linpackd.f"
      DO 8 I=1,II8
# 671 "linpackd.f"
            Y(I) = (((((((Y(I) + DD15 * M(I,II4-7)) + DD14 * M(I,II4-6))
     X              + DD13 * M(I,II4-5)) + DD12 * M(I,II4-4)) + DD11 * M
     X             (I,II4-3)) + DD10 * M(I,II4-2)) + DD9 * M(I,II4-1)) +
     X              DD8 * M(I,II4)
# 676 "linpackd.f"
    8       CONTINUE
# 670 "linpackd.f"
C$DOACROSS SHARE(II8,N1,Y,DD15,II4,M,DD14,DD13,DD12,DD11,DD10,DD9,DD8),
C$&  LOCAL(I)
# 670 "linpackd.f"
         DO 9 I=II8+1,N1,2
# 671 "linpackd.f"
            Y(I) = (((((((Y(I) + DD15 * M(I,II4-7)) + DD14 * M(I,II4-6))
     X              + DD13 * M(I,II4-5)) + DD12 * M(I,II4-4)) + DD11 * M
     X             (I,II4-3)) + DD10 * M(I,II4-2)) + DD9 * M(I,II4-1)) +
     X              DD8 * M(I,II4)
# 671 "linpackd.f"
            Y(I+1) = (((((((Y(I+1) + DD15 * M(I+1,II4-7)) + DD14 * M(I+1
     X             ,II4-6)) + DD13 * M(I+1,II4-5)) + DD12 * M(I+1,II4-4)
     X             ) + DD11 * M(I+1,II4-3)) + DD10 * M(I+1,II4-2)) + DD9
     X              * M(I+1,II4-1)) + DD8 * M(I+1,II4)
# 676 "linpackd.f"
    9    CONTINUE
# 677 "linpackd.f"
      endif
c
c   main loop - groups of sixteen vectors
c
      JMIN = MOD (N2, 16) + 16
# 683 "linpackd.f"
C$DOACROSS SHARE(N1,Y,JMIN,N2,X,M),LOCAL(DD16,I,J)
# 683 "linpackd.f"
         DO 11 I=1,N1
# 683 "linpackd.f"
         DD16 = Y(I)
# 682 "linpackd.f"
      DO 60 J=JMIN,N2,16
# 684 "linpackd.f"
            DD16 = (((((((((((((((DD16 + X(J-15) * M(I,J-15)) + X(J-14)
     X             * M(I,J-14)) + X(J-13) * M(I,J-13)) + X(J-12) * M(I,J
     X             -12)) + X(J-11) * M(I,J-11)) + X(J-10) * M(I,J-10)) +
     X              X(J-9) * M(I,J-9)) + X(J-8) * M(I,J-8)) + X(J-7) * M
     X             (I,J-7)) + X(J-6) * M(I,J-6)) + X(J-5) * M(I,J-5)) +
     X             X(J-4) * M(I,J-4)) + X(J-3) * M(I,J-3)) + X(J-2) * M(
     X             I,J-2)) + X(J-1) * M(I,J-1)) + X(J) * M(I,J)
# 694 "linpackd.f"
   60 continue
      Y(I) = DD16
# 693 "linpackd.f"
   11    CONTINUE
# 695 "linpackd.f"
      return
      end
 
 
